{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center><h1> Generalized linear models</h1></center>\n",
    "\n",
    "# 1. The exponential family\n",
    "We have now encountered a wide variery of probability distributions: the Gaussian, the Bernoulli, the Student t, the uniform, the gamma, etc. It turns out that most of these are members of a broader class of distributions known as the **exponential family**.\n",
    "\n",
    "We will see how we can easily use any member of the exponential family as a class-conditional density in order to make a generative classifier. In addition, we will discuss how to build discriminative models, where the response variable has an exponential family distribution,whose mean parameters is a linear function of the inputs, which is known as a generalized linear model, and generalizes the idea of logistic regression to other kind of response variables.\n",
    "\n",
    "## 1.1 Definition\n",
    "A pdf  or pmf $p(\\vec{x}|\\vec{\\theta})$ for $\\vec{x}=(x_1,x_2,\\ldots,x_m) \\in \\mathcal{X}^m$ and $\\vec{\\theta} \\in \\Theta \\subseteq R^d$ is said to be in the exponential family if it is of the form\n",
    "\\begin{align}\n",
    "p(\\vec{x}|\\vec{\\theta})&=\\frac{1}{Z(\\vec{\\theta})}h(\\vec{x})exp \\left[\\vec{\\theta}^T\\phi(\\vec{x}) \\right] \\\\\n",
    "                       &=h(\\vec{x})exp \\left[\\vec{\\theta}^T\\phi(\\vec{x})-A(\\vec{\\theta}) \\right]\n",
    "\\end{align}\n",
    "where\n",
    "\\begin{align}\n",
    "Z(\\vec{\\theta}) &= \\int h(\\vec{x})exp \\left[\\vec{\\theta}^T\\phi(\\vec{x}) \\right] d\\vec{x} \\\\\n",
    "A(\\vec{\\theta}) &=log\\,Z(\\vec{\\theta})\n",
    "\\end{align}\n",
    "Here $\\vec{\\theta}$ is called the **natural parameters** or **canonical parameters**,$\\phi(\\vec{x})$ is called a vector of **sufficient statistics**,$Z(\\vec{\\theta})$ is called the **partition function**,$A(\\vec{\\theta})$ is called the **cumulant function**, $h(\\vec{x})$ is the **scaling constant**,which is often 1\n",
    "\n",
    "## 1.2 Some examples\n",
    "### 1.2.1 Bernoulli distribution\n",
    "The Bernoulli for $x \\in \\lbrace 0,1 \\rbrace$ can be writen as\n",
    "\\begin{align}\n",
    "p(x|\\mu) &=\\mu^x(1-\\mu)^{1-x} \\\\\n",
    "         &=exp[x\\, log(\\mu)+(1-x)\\,log(1-\\mu)] \\\\\n",
    "         &=exp[x\\,log\\frac{\\mu}{1-\\mu}+log(1-\\mu)]\\\\\n",
    "         &=(1-\\mu) exp[log \\left(\\frac{\\mu}{1-\\mu} \\right) x] \\\\\n",
    "\\end{align}\n",
    "Now we have \n",
    "\\begin{align}\n",
    "\\phi(x) &=x  \\\\\n",
    "\\theta        &=log\\,\\left(\\frac{\\mu}{1-\\mu} \\right)\n",
    "\\end{align}\n",
    "which is the **sufficient statistics** and **natural parameters** for the Bernoulli distribution. We can recover the mean parameter $\\mu$ from the natural parameters using\n",
    "$$\n",
    "\\mu=sigm(\\theta)=\\frac{1}{1+e^{-\\theta}}\n",
    "$$\n",
    "### 1.2.2 Multinoulli\n",
    "The Multinoulli for $x \\in \\lbrace 1,\\ldots,C \\rbrace$ can be writen as \n",
    "\\begin{align}\n",
    "Cat(x|\\vec{\\mu})&=\\prod_{c=1}^C \\mu_c^{\\mathbb{1}(x=c)} \\\\\n",
    "              &=exp[\\sum_{c=1}^C \\mathbb{1}(x=c)\\,log\\,\\mu_c]   \\\\\n",
    "              &=exp[\\sum_{c=1}^{C-1} \\mathbb{1}(x=c) \\,log\\,\\mu_c +\\left(1-\\sum_{c=1}^{C-1} \\mathbb{1}(x=c) \\right) \\,log\\, \\mu_C] \\\\\n",
    "              &=exp[\\sum_{c=1}^{C-1} \\mathbb{1}(x=c) \\,log \\, \\frac{\\mu_c}{\\mu_C} + log\\,\\mu_C ]  \\\\\n",
    "              &=(1-\\sum_{c=1}^{C-1} \\mu_c)exp[\\sum_{c=1}^{C-1} \\mathbb{1}(x=c) \\,log \\, \\frac{\\mu_c}{\\mu_C}]\n",
    "\\end{align}\n",
    "Now we have \n",
    "\\begin{align}\n",
    "\\phi(x) &=\\lbrace \\mathbb{1}(x=1),\\ldots,\\mathbb{1}(x=C-1) \\rbrace \\\\\n",
    "\\vec{\\theta}  &=\\lbrace log\\,\\frac{\\mu_1}{\\mu_C},\\ldots,log\\,\\frac{\\mu_{C-1}}{\\mu_C} \\rbrace \\\\\n",
    "\\end{align}\n",
    "which is the **sufficient statistics** and **natural parameters** for the Multinoulli distribution. The form of the natural parameters for the multinoulli distribution is similar to the Bernoulli distribution.\n",
    "\n",
    "### 1.2.3 Univariate Gaussian \n",
    "The univariate Gaussian can be written in exponential family form as follows\n",
    "\\begin{align}\n",
    "\\mathcal{N}(x|\\mu,\\sigma^2) &= \\frac{1}{(2\\pi\\sigma^2)^{\\frac{1}{2}}}\n",
    "exp\\left[-\\frac{(x-\\mu)^2}{2\\sigma^2} \\right] \\\\\n",
    "                            &= \\frac{1}{(2\\pi\\sigma^2)^{\\frac{1}{2}}}exp \\left[-\\frac{x^2}{2\\sigma^2}+\\frac{\\mu x}{\\sigma^2}-\\frac{\\mu^2}{2\\sigma^2} \\right]\n",
    "\\end{align}\n",
    "Now we have \n",
    "\\begin{align}\n",
    "\\phi(x)  &=\\lbrace x^2,x \\rbrace \\\\\n",
    "\\vec{\\theta}  &=\\lbrace \\frac{-1}{2\\sigma^2},\\frac{\\mu}{\\sigma^2} \\rbrace \\\\\n",
    "\\end{align}\n",
    "which is the **sufficient statistics** and **natural parameters** for the Univariate Gaussian distribution.\n",
    "\n",
    "## 1.3 MLE for the exponential family\n",
    "The likelihood of an exponential family model has the form\n",
    "$$\n",
    "p(D|\\vec{\\theta})=\\left(\\frac{1}{Z(\\vec{\\theta})}\\right)^N \\left[\\prod_{i=1}^N h(\\vec{x}_i) \\right] exp\\left[ \\vec{\\theta}^T \\sum_{i=1}^N \\phi(\\vec{x}_i) \\right]\n",
    "$$\n",
    "We see that the sufficient statistics are $N$ and\n",
    "$$\n",
    "\\vec{\\phi}(D)=[\\sum_{i=1}^N \\phi_1(x_i),\\ldots,\\sum_{i=1}^N \\phi_K(x_i)]\n",
    "$$\n",
    "The log-likelihood is that\n",
    "$$\n",
    "log\\, p(D|\\vec{\\theta})=\\vec{\\theta}^T \\vec{\\phi}(D)-N\\,log\\,Z(\\vec{\\theta})-\\sum_{i=1}^N log\\,h(\\vec{x}_i)\n",
    "$$\n",
    "The derivative of the log-likelihood is that\n",
    "\\begin{align}\n",
    "\\frac{\\partial}{\\partial \\vec{\\theta}} log\\, p(D|\\vec{\\theta}) &=\\vec{\\phi}(D)-N\\,\\frac{\\partial}{\\partial \\vec{\\theta}} log\\,Z(\\vec{\\theta}) \\\\\n",
    "&=\\vec{\\phi}(D)-N\\,\\frac{1}{Z(\\vec{\\theta})}\\frac{\\partial}{\\partial \\vec{\\theta}}Z(\\vec{\\theta}) \\\\\n",
    "&=\\vec{\\phi}(D)-N\\,\\frac{1}{Z(\\vec{\\theta})}\\frac{\\partial}{\\partial \\vec{\\theta}} \\int h(\\vec{x})exp\\{\\vec{\\theta}^T\\phi(\\vec{x})\\}d\\vec{x} \\\\\n",
    "&=\\vec{\\phi}(D)-N\\,\\frac{1}{Z(\\vec{\\theta})}\\int \\frac{\\partial}{\\partial \\vec{\\theta}} h(\\vec{x})exp\\{\\vec{\\theta}^T\\phi(\\vec{x})\\}d\\vec{x} \\\\\n",
    "&=\\vec{\\phi}(D)-N\\,\\frac{1}{Z(\\vec{\\theta})}\\int \\phi(\\vec{x})  h(\\vec{x})exp\\{\\vec{\\theta}^T\\phi(\\vec{x})\\}d\\vec{x} \\\\\n",
    "&=\\vec{\\phi}(D)-N\\,\\int \\phi(\\vec{x}) \\frac{1}{Z(\\vec{\\theta})}  h(\\vec{x})exp\\{\\vec{\\theta}^T\\phi(\\vec{x})\\}d\\vec{x} \\\\\n",
    "&=\\vec{\\phi}(D)-N\\,\\mathbb{E}[\\phi(\\vec{x})]\n",
    "\\end{align}\n",
    "So the MLE means the empirical average of the sufficient statistics must equal the model’s theoretical expected sufficient statistics, which means\n",
    "$$\n",
    "\\mathbb{E}[\\phi(\\vec{x})]=\\frac{1}{N}\\sum_{i=1}^N\\phi(\\vec{x}_i)\n",
    "$$\n",
    "Think about this, which is just a **moment matching** process.\n",
    "## 1.4 Bayes for the exponential family\n",
    "The likelihood of the exponential family is given by\n",
    "$$\n",
    "p(D|\\vec{\\theta}) \\propto g(\\vec{\\theta})^N exp(\\vec{\\theta}^T\\vec{s}_N)\n",
    "$$\n",
    "where \\begin{align}\n",
    "g(\\vec{\\theta}) &=\\frac{1}{Z(\\vec{\\theta})}\\\\\n",
    "\\vec{s}_N  &=\\sum_{i=1}^N \\phi(\\vec{x}_i)\n",
    "\\end{align}\n",
    "The natural conjugate prior has the form\n",
    "\\begin{align}\n",
    "p(\\vec{\\theta}|\\nu_0,\\vec{\\tau}_0) &\\propto g(\\vec{\\theta})^{\\nu_0} exp \\left(\\vec{\\theta}^T \\vec{\\tau}_0 \\right)\\\\\n",
    "                                   &=g(\\vec{\\theta})^{\\nu_0} exp \\left(\\nu_0 \\vec{\\theta}^T \\vec{\\bar{\\tau}}_0 \\right) \\\\\n",
    "                                   &=p(\\vec{\\theta}|\\nu_0,\\vec{\\bar{\\tau}})\n",
    "\\end{align}\n",
    "The $\\nu_0$ means the size of  prior pseudo-data and $\\vec{\\bar{\\tau}}_0$ means the mean of the sufficient statistics on this pseudo-data.The posterior is given by\n",
    "\\begin{align}\n",
    "p(\\vec{\\theta}|D) &\\propto p(D|\\vec{\\theta})p(\\vec{\\theta}|\\nu_0,\\vec{\\tau}_0) \\\\\n",
    "                  &=g(\\vec{\\theta})^{\\nu_0+N} exp \\left(\\vec{\\theta}^T(\\nu_0 \\vec{\\bar{\\tau}}_0+\\vec{s}_N) \\right)\n",
    "\\end{align}\n",
    "\n",
    "## 1.5 Maximum entropy derivation of the exponential family\n",
    "The exponential family distribution is that makes the least number of assumptions about the data ,subject to a specific set of user-specified constraints. Suppose all we known is the expected values of certain features or functions:\n",
    "$$\n",
    "\\sum_x f_k(x)p(x)=F_k,\\quad k=1,2,\\ldots,K\n",
    "$$\n",
    "where $F_k$ are known constants, and $f_k(x)$ si an arbitrary function.The principle of **maximum\n",
    "entropy** or **maxent** says we should pick the distribution with maximum entropy (closest to\n",
    "uniform), subject to the constraints that the moments of the distribution match the empirical\n",
    "moments of the specified functions. To maximum entropy subject to the constraints, we need to use Lagrange multipliers as follows\n",
    "$$\n",
    "\\mathcal{J}(p,\\lambda)=-\\sum_xp(x)log\\,p(x)+\\lambda_0\\left(1-\\sum_xp(x)\\right)+\\sum_k \\lambda_k\\left(F_k-\\sum_x f_k(x)p(x)\\right)\n",
    "$$\n",
    "We can use the calculus of variations to take derivatives wrt the distribution $p(x)$, then we have\n",
    "$$\n",
    "\\frac{\\partial \\mathcal{J}}{\\partial p(x)}=-1-log\\,p(x)-\\lambda_0-\\sum_k \\lambda_k f_k(x)\n",
    "$$\n",
    "Setting $\\frac{\\partial \\mathcal{J}}{\\partial p(x)}=0$ yields\n",
    "$$\n",
    "p(x)=\\frac{1}{Z} exp\\left(-\\sum_k \\lambda_k f_k(x)\\right)\n",
    "$$\n",
    "which is the exponential family distribution, also known as the **Gibbs distribution**.\n",
    "\n",
    "# 2. Generalized linear models(GLMs)\n",
    "## 2.1 Basic\n",
    "These are models in which the output density is in the exponential family, and in which the mean parameters are a linear combination of the inputs, passed through a possibly nonlinear function.\n",
    "\n",
    "The exponential family is written as \n",
    "\\begin{align}\n",
    "p(\\vec{y}|\\vec{\\theta}) &=h(\\vec{y})exp \\left[\\vec{\\theta}^T\\vec{y}-A(\\vec{\\theta}) \\right] \\\\\n",
    "                        &=exp \\left[\\vec{\\theta}^T\\vec{y}-A(\\vec{\\theta})+C(\\vec{y}) \\right]\n",
    "\\end{align}\n",
    "where $C(\\vec{y})=log\\,h(\\vec{y})$\n",
    "\n",
    "More generally, one sometimes introduces an extra parameter $\\sigma^2$, called the dispersion parameter, to control the shape of $p(\\vec{y}|\\vec{\\theta})$ as follows\n",
    "$$\n",
    "p(\\vec{y}|\\vec{\\theta},\\sigma^2)=exp \\left[\\frac{\\vec{\\theta}^T\\vec{y}-A(\\vec{\\theta})}{\\sigma^2}+C(\\vec{y},\\sigma^2) \\right]\n",
    "$$\n",
    "where $\\sigma^2$ is the **dispersion parameter** and $\\vec{\\theta}$ is the **natural parameter**, $A$ is the partition function. \n",
    "\n",
    "We now define a linear function of the inputs $\\vec{x}$\n",
    "$$\n",
    "\\vec{\\eta}=\\mathbf{W}\\vec{x}\n",
    "$$\n",
    "Then,we make the mean of the distribution be some invertible monotonic function of $\\vec{\\eta}$, which known as the **mean function**, is denoted by $g^{-1}$\n",
    "$$\n",
    "\\vec{\\mu}=\\mathbb{E}[\\vec{y}|\\vec{\\theta},\\sigma^2]=g^{-1}(\\vec{\\eta})=g^{-1}(\\mathbf{W}\\vec{x})\n",
    "$$\n",
    "The inverse of the mean function ,namely $g$ ,is called the **link function**.\n",
    "\n",
    "One particularly simple form of link function is to make $\\vec{\\theta}=\\vec{\\eta}$,which is called **canonical link function**.\n",
    "### 2.1.1 Linear regression\n",
    "In the linear regression case, we care much more attention on the mean $\\mu$ than the variance $\\sigma^2$. We can view $\\sigma^2$ as a dispersion parameter.\n",
    "\\begin{align}\n",
    "p(y|\\mu,\\sigma^2) &=\\mathcal{N}(y|\\mu,\\sigma^2) \\\\\n",
    "                  &=exp \\left[-\\frac{(y-\\mu)^2}{2\\sigma^2}-\\frac{1}{2}log\\,\\left(2\\pi\\sigma^2 \\right) \\right] \\\\\n",
    "                  &=exp \\left[\\frac{y\\mu-\\frac{\\mu^2}{2}}{\\sigma^2}-\\frac{1}{2}\\left(\\frac{y^2}{\\sigma^2}+log\\,\\left(2\\pi\\sigma^2 \\right) \\right) \\right]\n",
    "\\end{align}\n",
    "The **natural parameter** is $\\theta=\\mu=\\vec{w}^T\\vec{x}$, which is the ordinary linear regression case.\n",
    "\n",
    "### 2.1.2 Binomial regression\n",
    "We have \n",
    "\\begin{align}\n",
    "p(y|\\mu) &=Bin(y|N,\\mu) \\\\\n",
    "         &=\\binom{N}{y}\\mu^y(1-\\mu)^{N-y} \\\\\n",
    "         &=exp \\left[log\\,\\binom{N}{y}+ylog\\,\\mu+(N-y)log\\,(1-\\mu) \\right] \\\\\n",
    "         &=exp \\left[ylog\\,\\frac{\\mu}{1-\\mu}+Nlog\\,(1-\\mu)+log\\,\\binom{N}{y} \\right] \\\\\n",
    "\\end{align}\n",
    "The **natural parameter** is \n",
    "$$\n",
    "\\theta=log\\,\\frac{\\mu}{1-\\mu}=\\vec{w}^T\\vec{x}\n",
    "$$\n",
    "### 2.1.3 Binary classification\n",
    "In (binary) logistic regression, we use a model of the form\n",
    "$$\n",
    "p(y=1|\\vec{x},\\vec{w})=sigm(\\vec{w}^T\\vec{x})\n",
    "$$\n",
    "In general, we can write\n",
    "$$\n",
    "p(y=1|\\vec{x},\\vec{w})=g^{-1}(\\vec{w}^T\\vec{x})\n",
    "$$\n",
    "for any function $g^{-1}$ which maps $[-\\infty,\\infty]$ to $[0,1]$. Several possible mean functions are listed below.\n",
    "\n",
    "| Name     |          Formula &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;&nbsp; &nbsp; &nbsp; &nbsp; &nbsp;&nbsp; &nbsp; &nbsp; &nbsp; &nbsp;\n",
    "|----------|:----------------------------------------:|\n",
    "| Logistic |         $g^{-1}(\\eta)=sigm(\\eta)$        |\n",
    "| Probit   |         $g^{-1}(\\eta)=\\Phi(\\eta)$        |\n",
    "| Log-log  |$g^{-1}(\\eta)=exp\\left(-exp(-\\eta)\\right)$|\n",
    "## 2.2 ML and MAP estimation\n",
    "The likelihood has the following form\n",
    "\\begin{align}\n",
    "p(D|\\vec{w}) &=\\prod_{i=1}^N exp \\left[\\frac{\\theta_i y_i-A(\\theta_i)}{\\sigma^2}+C(y_i,\\sigma^2) \\right] \\\\\n",
    "\\theta_i     &=\\vec{w}^T\\vec{x}_i\n",
    "\\end{align}\n",
    "We can compute the gradient respect to $\\vec{w}$ as follows\n",
    "\\begin{align}\n",
    "\\frac{\\partial log\\, p(D|\\vec{w})}{\\partial \\vec{w}} &=\\sum_{i=1}^N\\frac{\\partial log\\, p(y_i|\\vec{w},\\vec{x_i})}{\\partial \\theta_i}\\frac{\\partial \\theta_i}{ \\partial \\vec{w}} \\\\\n",
    "&=\\sum_{i=1}^N \\frac{y_i-A'(\\theta_i)}{\\sigma^2}\\frac{\\partial \\theta_i}{ \\partial \\vec{w}} \\\\\n",
    "&=\\sum_{i=1}^N \\frac{y_i-A'(\\theta_i)}{\\sigma^2}\\vec{x}_i\n",
    "\\end{align}\n",
    "in which $y_i-A'(\\theta_i)$ takes the error form. The gradient form is similar to the logistic regression. It is straightforward to modify the above procedure to perform MAP estimation with a Gaussian prior"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [default]",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
